Problem Set 1 Key: Statistics Fundamentals

PLS 206 — Applied Multivariate Modeling

Author

Grey Monroe

Published

January 1, 2026

Tip

Try the problems yourself before consulting this key. Working through errors is where most of the learning happens.

Part 1: Create and Manipulate a Data Frame

Q1. set.seed() fixes the starting point of R’s random number generator. Using the same seed ensures that code producing random values gives identical results every time it is run — essential for reproducible analyses.

?set.seed
# or
help(set.seed)

Q2. Create df and find the sum of all values.

set.seed(1522)
df <- data.frame(
  yield = rnorm(100, mean = 100, sd = 5),
  temp  = rnorm(100, mean = 10,  sd = 1)
)
sum(df)

Q3. Mean of yield:

mean(df$yield)

Q4. 8th value of yield:

df$yield[8]

Part 2: Import and Examine the Wine Dataset

wine <- read.csv("winedata.csv", header = TRUE)
str(wine)

Q5. Dimensions:

dim(wine)   # 178 rows, 14 columns

Q6. Object type:

class(wine)   # "data.frame"

Q7. Variable types:

str(wine)
# a) Categorical:          Cultivar
# b) Quantitative integer: Mg or Proline
# c) Quantitative numeric: Alcohol, Ash, Color, Flav, OD, etc.

Part 3: Identify and Fix Common Coding Errors

Q8. Bug: wrong capitalization (DF vs df).

# Buggy
DF[1:5, 1:2]

# Fixed
df[1:5, 1:2]

Q9. Bug: missing comma between arguments in rbind().

# Buggy
dataframe <- rbind(vector.a vector.b)

# Fixed
vector.a  <- c(0.2, 0.5, 0.8, 0.01, 0.03)
vector.b  <- c(12, 7, 5, 4, 14)
dataframe <- rbind(vector.a, vector.b)

Q10. Bug: missing $ operator and incorrect bracket syntax.

# Buggy
color5 <- wine[wine Color >= 5]

# Fixed
color5 <- subset(wine, Color >= 5)
# or equivalently:
color5 <- wine[wine$Color >= 5, ]

Part 4: Calculate Basic Summary Statistics

Q11. Mean of all quantitative columns — Proline has the highest mean.

wine_means <- apply(wine[, 2:14], 2, mean)
wine_means
which.max(wine_means)   # Proline

Q12. SD of all quantitative columns — NonFlavPhenols has the lowest SD.

wine_sds <- apply(wine[, 2:14], 2, sd)
wine_sds
which.min(wine_sds)   # NonFlavPhenols

Part 5: Subset the Dataset

Q13. Reduce to selected columns and display the top rows.

# By column name (preferred — more readable)
wine_reduce <- wine[, c("Cultivar", "Alcohol", "Ash", "Color", "Flav", "Mg", "OD")]

# Equivalent using column indices
wine_reduce <- wine[, c(1, 2, 4, 5, 6, 9, 11)]

head(wine_reduce)

Part 6: Examine Correlations Between Variables

Q14 & Q15. Scatter plot matrix — strongest correlations (by absolute value):

  1. Flav & OD: r = 0.787
  2. Color & Alcohol: r = 0.546
  3. OD & Color: r = −0.429
library(ggplot2)
library(GGally)

ggpairs(wine_reduce, columns = 2:7)

Part 7: Examine Data by Cultivar

Q16. Sample size per cultivar:

table(wine_reduce$Cultivar)
# Barbera: 48   Barolo: 59   Gringnolino: 71

Q17. Colored scatter plot matrix by cultivar:

ggpairs(wine_reduce, aes(col = Cultivar))

Q18. Mean by cultivar:

aggregate(. ~ Cultivar, data = wine_reduce, FUN = mean)

# a) Clearly different across cultivars: Alcohol, Color, Flav, OD
# b) Mostly overlapping:                Ash, Mg

Part 8: Flav vs. OD

Q19. The overall correlation between Flav and OD is positive (~0.79), but this masks very different within-cultivar patterns:

# Overall
cor(wine_reduce$Flav, wine_reduce$OD)

# By cultivar
by(wine_reduce[, c("Flav", "OD")], wine_reduce$Cultivar,
   function(x) cor(x$Flav, x$OD))

# Barbera:     r ≈ −0.430
# Barolo:      r ≈ −0.089
# Gringnolino: r ≈ +0.580

The overall positive correlation is driven primarily by Gringnolino. Within Barbera and Barolo the relationship is flat or slightly negative — a clear example of Simpson’s Paradox, where a group-level trend reverses or disappears within subgroups.

Part 9: Uncertainty and Comparing Two Groups

Q20. Barolo alcohol: n = 59, mean = 13.74, SD = 0.462, SE = 0.0602, 95% CI = 13.62 to 13.87. t.test() gives the same interval.

x     <- wine$Alcohol[wine$Cultivar == "barolo"]
n     <- length(x)
xbar  <- mean(x)
s     <- sd(x)
se    <- s / sqrt(n)
tcrit <- qt(0.975, df = n - 1)

c(n = n, mean = xbar, sd = s, se = se,
  lower = xbar - tcrit * se, upper = xbar + tcrit * se)

t.test(x)$conf.int

Q21.

    1. The SD (0.46 % alcohol) is the spread of individual wines around the mean. The SE (0.06) is the spread we’d expect in the sample mean if we repeated the study, so it describes the precision of our estimate.
    1. The SD should stay about the same, because it’s a property of the wines, not of the sample size. The SE should roughly halve, since SE = s/√n and √4 = 2.

Q22. Barolo mean Proa = 1.90, Gringnolino = 1.63. Difference = 0.27 (95% CI 0.09 to 0.45), Welch t = 3.01, df ≈ 123.7, p = 0.0032.

Example interpretation: Barolo wines contained on average 0.27 units more proanthocyanins than Gringnolino wines (95% CI 0.09–0.45; Welch t = 3.01, p = 0.003).

two <- subset(wine, Cultivar %in% c("barolo", "gringnolino"))

ggplot(two, aes(Cultivar, Proa)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.1, alpha = 0.6) +
  theme_classic()

t.test(Proa ~ Cultivar, data = two)

Q23. The permutation p-value should be close to the t-test’s (≈ 0.003–0.004; with set.seed(206) it is 0.0036). Exact values vary with the seed. Both methods ask how often chance alone would produce a difference this large.

set.seed(206)
observed <- mean(two$Proa[two$Cultivar == "barolo"]) -
            mean(two$Proa[two$Cultivar == "gringnolino"])

null_diffs <- replicate(10000, {
  shuffled <- sample(two$Cultivar)
  mean(two$Proa[shuffled == "barolo"]) - mean(two$Proa[shuffled == "gringnolino"])
})

mean(abs(null_diffs) >= abs(observed))

Part 10: More Than Two Groups and Effect Sizes

Q24. F = 233.9 on 2 and 175 df, p < 2e-16. Tukey: all three pairs differ (Barolo > Gringnolino > Barbera; adjusted p < 0.001 for each).

a_flav <- aov(Flav ~ Cultivar, data = wine)
summary(a_flav)
TukeyHSD(a_flav)

Q25. The residuals depart from normality (Shapiro–Wilk W = 0.933, p = 2.3e-07; the Q–Q plot bends at the tails). The Kruskal–Wallis test is also highly significant (χ² = 130.5, df = 2, p < 2.2e-16), so the conclusion doesn’t change. With a difference this large and ~60 wines per group, the ANOVA is robust to the violation.

qqnorm(residuals(a_flav)); qqline(residuals(a_flav))
shapiro.test(residuals(a_flav))
kruskal.test(Flav ~ Cultivar, data = wine)

Q26.

    1. Cohen’s d, Barolo vs. Gringnolino: Proa d ≈ 0.51 (medium); Alcohol d ≈ 2.90 (enormous).
    1. Eta-squared: Flav η² ≈ 0.73 (cultivar explains ~73% of the variance); Mg η² ≈ 0.12 (~12%).
    1. A p-value reflects both effect size and sample size, so with enough data a small difference can be highly significant. The Mg ANOVA (p ≈ 9e-06) and the Flav ANOVA (p < 2e-16) are both “significant”, but cultivar explains six times more of the variance in Flav. Effect sizes measure how big the difference is, which is the biologically meaningful question.
cohens_d <- function(x1, x2) {
  n1 <- length(x1); n2 <- length(x2)
  pooled_sd <- sqrt(((n1 - 1) * var(x1) + (n2 - 1) * var(x2)) / (n1 + n2 - 2))
  abs(mean(x1) - mean(x2)) / pooled_sd
}
b <- two$Cultivar == "barolo"
cohens_d(two$Proa[b], two$Proa[!b])
cohens_d(two$Alcohol[b], two$Alcohol[!b])

eta_sq <- function(model) {
  ss <- summary(model)[[1]][, "Sum Sq"]
  ss[1] / sum(ss)
}
eta_sq(a_flav)
eta_sq(aov(Mg ~ Cultivar, data = wine))